/*
 * FVoronoi.cpp
 *
 *  Created on: May 13, 2012
 *      Author: marchi
 */

#include "FVoronoi.h"


const int NTYPES=3;
struct VolAreaInd{
	float Vol,area[NTYPES];
	int n;
	int Type;
	bool bContact;
	VolAreaInd():Vol(0),n(-1),Type(-1),bContact(false){for(int o=0;o<NTYPES;o++) area[o]=0.0;};
};
class ContactSolute{
public:
	bool operator()(VolAreaInd x){return (!x.bContact);}
};
class Solute{
public:
	bool operator()(VolAreaInd x){return x.Type != Slt;}
};
class Solvent{
public:
	bool operator()(VolAreaInd x){return x.Type != Sol;}
};

int FVoronoi::nstat=0;
double FVoronoi::dx=1.0;
double FVoronoi::dvx=0.07;
double FVoronoi::fVol=1000.0;
double FVoronoi::fVol2=fVol*fVol;
int FVoronoi::nx=static_cast<int> (300.0/dx);
int FVoronoi::nvx=static_cast<int> (20.0/dvx);

void FVoronoi::Stats(double frame){
	if(!nstat){
		double cut=(co[XX][XX]+co[YY][YY]+co[ZZ][ZZ])/6.0;
		nvx=static_cast<int> (cut*10.0/dvx);
		hist0=vector<double>(nx,0.0);
		vv0R=vector<double>(nvx,0.0);
		gofr=vector<double>(nvx,0.0);
		vR=vector<double>(nvx,0.0);
		nvv0R=vector<int>(nvx,0);
	}
	nstat++;
	VolumeCell+=(co[XX][XX]*(co[YY][YY]*co[ZZ][ZZ]-co[ZZ][YY]*co[YY][ZZ])
			  -co[YY][XX]*(co[XX][YY]*co[ZZ][ZZ]-co[ZZ][YY]*co[XX][ZZ])
			  +co[ZZ][XX]*(co[XX][YY]*co[YY][ZZ]-co[YY][YY]*co[XX][ZZ]));
	double fact=1.0/static_cast<double> (nresid);
	for(int n=0;n<nresid;n++){
		int mx=static_cast<int> (1000.0*Vols[n]/dx);
		if(mx < nx) hist0[mx]+=1.0;
	}
	int o=0;
	for(int n=0;n<nresid;n++){
		string labs=ResidueCM::getlabels(n);
		if(labs.find("SOL") == string::npos) continue;
		v0+=fact*Vols[n]*fVol;
		for(int m=n;m<nresid;m++){
			string labs1=ResidueCM::getlabels(m);
			if(labs1.find("SOL") == string::npos) continue;
			double rsp=10.0*Cxcm[n].Dist(Cxcm[m],co,oc);
/*
			if(rsp < 3 && rsp != 0.0){
				cout << n << " " << m << " " << rsp <<endl;
				cout << Cxcm[n][XX] << " " << Cxcm[n][YY] << " " << Cxcm[n][ZZ] << " " << endl;
				cout << Cxcm[m][XX] << " " << Cxcm[m][YY] << " " << Cxcm[m][ZZ] << " " << endl;
				exit(1);
			}
*/

			int myx=static_cast<int> (rsp/dvx);
			if(myx < nvx) {
				vv0R[myx]+=Vols[n]*Vols[m]*fVol2;
				vR[myx]+=Vols[m]*fVol;
				nvv0R[myx]++;
			}
		}
	}
}
void FVoronoi::WriteStats(Parallel::MyMPI * comm){
	v0/=static_cast<double>(nstat);
	VolumeCell/=static_cast<double>(nstat);
	for(int n=0;n<nvx;n++){
		gofr[n]=static_cast<double> (nvv0R[n])/static_cast<double>(nstat*nresid);
		double R=dvx*static_cast<double>(n);
		double dvol=2.0*PI*R*R*dvx;
		if(dvol) gofr[n]/=dvol/v0;
		else gofr[n]=0.0;
	}

	double Norm;
	vector<double> gofr0(nvx),delta2(nvx);
	for(int n=0;n<nvx;n++){
		if(nvv0R[n]){
			vR[n]/=static_cast<double>(nvv0R[n]);
			vv0R[n]/=static_cast<double>(nvv0R[n]);
		} else{
			vR[n]=0.0;
			vv0R[n]=0.0;
		}
		if(n==0) Norm=vv0R[n]-v0*vR[n];
		double avg=(vv0R[n]/(v0*vR[n])-1.0);
		delta2[n]=avg;
		gofr0[n]=gofr[n];
	}
	comm->ReduceSum(delta2);
	comm->ReduceSum(gofr0);
	double nMPI=static_cast<double> (comm->Get_Size());
	for(int n=0;n<nvx;n++){
		delta2[n]/=nMPI;
		gofr0[n]/=nMPI;
		cout << dvx*static_cast<double>(n) << " " << delta2[n]*gofr0[n] << " " << delta2[n]<< " " << gofr0[n] <<  " " << endl;
	}
	double sum=0.0;
	for(int n=0;n<nx;n++){
			hist0[n]/=static_cast<double>(nstat);
		sum+=hist0[n];
	}

	sum*=dx;
	for(int n=0;n<nx;n++){
		hist0[n]/=sum;
		if(hist0[n]){
			cout << dx*static_cast<double>(n) << " " << hist0[n] << endl;
		}
	}

}
std::ofstream & operator<<(std::ofstream & fout , FVoronoi  & y){
	vector<VolAreaInd> tmp(y.nr);
	for(int o=0;o<y.nresid;o++) {
		vector<int> & cindex=ResidueCM::getind(o);
		bool bC=false;
		for(unsigned int ia=0;ia<cindex.size();ia++){
			int i=cindex[ia];
			tmp[i].Vol=static_cast<float> (y.Vol[i]*1000.0);
			for(unsigned int p=0;p<y.Neighs[i].size();p++){
				tmp[i].area[y.getTypes(y.Neighs[i][p])]+=static_cast<float> (y.Surface[i][p]*100.0);
			}
			if(tmp[i].area[0]) bC=true;
			tmp[i].n=i+1;
			tmp[i].Type=y.getTypes(i);
		}
		for(unsigned int ia=0;ia<cindex.size();ia++)
			tmp[cindex[ia]].bContact=bC;
	}
	vector<VolAreaInd> tmp0(tmp);
	tmp0.erase(std::remove_if(tmp0.begin(),tmp0.end(),Solute()),tmp0.end());
	tmp.erase(std::remove_if(tmp.begin(),tmp.end(),Solvent()),tmp.end());
	tmp.erase(std::remove_if(tmp.begin(),tmp.end(),ContactSolute()),tmp.end());
	vector<VolAreaInd> data;

	copy(tmp0.begin(),tmp0.end(),back_inserter(data));
	copy(tmp.begin(),tmp.end(),back_inserter(data));
	tmp.clear();
	tmp0.clear();

	int F4bytes=2*INT;
	int nstep=static_cast<int> (y.time);
	int nat=static_cast<int> (data.size());
	fout.write(reinterpret_cast<char *> (&F4bytes),sizeof(F4bytes));
		fout.write(reinterpret_cast<char *> (&nstep),sizeof(nstep));
		fout.write(reinterpret_cast<char *> (&nat),sizeof(nat));
	fout.write(reinterpret_cast<char *> (&F4bytes),sizeof(F4bytes));

	F4bytes=9*DOUBLE;
	double co[DIM][DIM];
	for(int o=0;o<DIM;o++)
		for(int p=0;p<DIM;p++)
			co[o][p]=y.co[o][p];

	fout.write(reinterpret_cast<char *> (&F4bytes),sizeof(F4bytes));
		fout.write(reinterpret_cast<char *> (co),sizeof(co));
	fout.write(reinterpret_cast<char *> (&F4bytes),sizeof(F4bytes));

	F4bytes=(INT+3*FLOAT)*nat;
	fout.write(reinterpret_cast<char *> (&F4bytes),sizeof(F4bytes));
	for(int n=0;n<nat;n++){
		int iat=data[n].n;
		float mydata[3]={data[n].Vol, data[n].area[0],data[n].area[1]};
		fout.write(reinterpret_cast<char *> (&iat),sizeof(iat));
		fout.write(reinterpret_cast<char *> (mydata),sizeof(mydata));
	}
	fout.write(reinterpret_cast<char *> (&F4bytes),sizeof(F4bytes));
	return fout;
}
FVoronoi::~FVoronoi() {
	// TODO Auto-generated destructor stub
}

